#!/usr/bin/env python

from osgeo import osr
from osgeo import gdal
from osgeo import ogr
from PIL import Image
from PIL import ImageDraw
from PIL import ImageFont
import base64
import StringIO
import math
import struct

reduction_tolerances = {
1:78271.52,
2:39135.76,
3:19567.88,
4:9783.94,
5:4891.97,
6:2445.98,
7:1222.99,
8:611.50,
9:305.75,
10:152.87,
11:76.44,
12:38.22,
13:19.11,
14:9.55,
15:4.78,
16:2.39,
17:1.19,
18:0.60,
19:0.30,
}

class georef:

  def __init__(self, upper_left, lower_right, lower_left, upper_right, width, height):
    self.upper_left = upper_left
    self.lower_right = lower_right
    self.upper_right = upper_right
    self.lower_left = lower_left
    merc_proj = georef.get_mercator_projection()
    def_proj = georef.get_default_projection()
    self.def_transform = osr.CoordinateTransformation(def_proj, merc_proj)
    geo = georef.get_geo_transform(upper_left, lower_right, width, height, self.def_transform)
    self.inv_geo = gdal.InvGeoTransform(geo)[1]
    self.img = Image.new("RGBA",(width,height),(0,0,0,0))
    self.width = width
    self.height = height
    self.merc_proj = merc_proj

  def get_bounding_box(self):
    return [self.lower_left[1],self.lower_left[0],self.upper_right[1],self.upper_right[0]]

  @staticmethod
  def make_georef(form):
    ll_lat = float(form["ll_lat"].value)
    ll_lon = float(form["ll_lon"].value)
    ur_lat = float(form["ur_lat"].value)
    ur_lon = float(form["ur_lon"].value)
    width = int(form["width"].value)
    height = int(float(form["height"].value))
    zoom = float(form["zoom"].value)
    ll = (ll_lat,ll_lon)
    ur = (ur_lat,ur_lon)
    ul = (ur[0],ll[1])
    lr = (ll[0],ur[1])
    gref = georef(ul,lr,ll,ur,width,height)
    return gref

  def get_image(self):
    f = StringIO.StringIO()
    self.img.save(f,"PNG")
    return f.getvalue()
    
  def get_base64(self):
    f = StringIO.StringIO()
    self.img.save(f,"PNG")
    return base64.b64encode(f.getvalue())

  def draw_polygon_from_geom(self,geom,color):
    polys = []
    for i in range(geom.GetGeometryCount()):
      poly = geom.GetGeometryRef(i).GetGeometryRef(0)
      poly_pts = map(lambda j: poly.GetPoint(j), range(poly.GetPointCount()))
      poly = self.draw_polygon(poly_pts,color)
      polys.append(poly)
    return polys

  def draw_polygon(self, points, color):
    dr = ImageDraw.Draw(self.img)
    points = map(lambda x: georef.get_pixel(x[0],x[1],self.inv_geo), points[:-1])
    dr.polygon(points, fill = color, outline = (0,0,0,255))
    del dr
    return points

  def draw_line(self, points, color):
    dr = ImageDraw.Draw(self.img)
    px_points = map(lambda pt: georef.get_pixel_from_lat_lon(pt[0],pt[1],self.inv_geo), points)
    unique_points = []
    [unique_points.append(i) for i in px_points if not unique_points.count(i)]
    dr.line(unique_points, fill = color, width = 3)
    dr.line(points, fill = color)
    del dr
    return unique_points

  def write_image(self, fname):
    self.img.save(fname)

  @staticmethod
  def get_mercator_projection():
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(900913)
    return proj

  @staticmethod
  def get_default_projection():
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(4326)
    return proj

  @staticmethod
  def get_geo_transform(upper_left, lower_right, width, height, transform):
    (ul_px, ul_py, z) = transform.TransformPoint(upper_left[1],upper_left[0])
    (lr_px, lr_py, z) = transform.TransformPoint(lower_right[1],lower_right[0])
    dx = abs(ul_px - lr_px)
    dy = abs(ul_py - lr_py)
    pixel_size_x = dx / float(width)
    pixel_size_y = dy / float(height)
    if ul_px < 0 and pixel_size_x < 0: # does this make sense?
      pixel_size_x *= -1             #this might only work for North America...
    if ul_py > 0 and pixel_size_y > 0:
      pixel_size_y *= -1
    return [ul_px,
              pixel_size_x,
              0,
              ul_py, #these values are in a weird order but this is correct
              0,
              pixel_size_y]

  @staticmethod
  def get_pixel_from_lat_lon(lat, lon, inv_geo, transform):
    (gx,gy,gz) = transform.TransformPoint(lon,lat)
    return georef.get_pixel(gx,gy,inv_geo)

  @staticmethod
  def get_pixel(x, y, inv_geo):
    (gx,gy) = gdal.ApplyGeoTransform(inv_geo, x, y)
    px = int(gx)
    py = int(gy)
    return (px,py)
